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Abstract 

A  numerical  method  is  developed  for  computing  steady  super¬ 
critical  flow  about  an  ellipse  at  zero  angle  of  attack.  The  flow 
is  assumed  to  be  two-dimensional,  inviscid,  isentropic,  and 
irrotational .  The  free  stream  Mach  number  lies  in  the  high 
subsonic  range  so  that  a  shock  wave  occurs  locally  near  the  body. 

The  full  potential  equations  are  solved  by  Telenin's  Method 
and  the  Method  of  Lines.  Smooth  interpolating  functions  are 
assumed  for  the  unknown  flow  variables  in  selected  coordinate 
directions.  The  resulting  set  of  ordinary  differential  equations 
is  then  integrated  away  from  or  along  the  body  depending  upon 
whether  the  flow  is  smooth  or  discontinuous.  Jump  conditions  of 
the  governing  equations  are  applied  across  the  shock  wave  so  that 
it  is  perfectly  sharp.  A  doublet  solution  for  flow  past  a  closed 
body  is  used  as  the  far  field  boundary  condition. 

Supercritical  flow  calculations  have  been  performed  for  ellipses 
with  thickness  ratio  of  0.2  and  0.4  at  various  free- stream  Mach 
numbers.  The  present  results  are  compared  with  the  shock-capturing 
method,  and  good  agreement  is  obtained. 
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1.0  Introduction 


Within  the  past  decade  there  has  been  strong  interest  in  transonic 
flow  research.  The  fact  that  commercial  jets  often  fly  at  transonic 
speeds  makes  it  very  desirable  to  have  methods  which  can  predict  air¬ 
foil  lift  and  drag  in  this  flow  regime. 

The  flow  is  called  transonic  if  both  subsonic  and  supersonic 
regions  are  present  in  the  field.  Although  most  airplanes  fly  at 
subsonic  speeds,  the  local  flow  velocities  often  become  supersonic 
at  the  top  of  the  wing.  In  a  typical  transonic  flow  field,  the 
embedded  supersonic  region  is  usually  terminated  by  means  of  a 
shock  wave. 

The  main  difficulties  in  transonic  flow  calculations  are  due  to 
the  inherent  nonlinearities  of  the  equations  governing  transonic  flow, 
and  the  fact  that  the  equations  change  type  within  the  solution 
domain,  from  elliptic  in  the  subsonic  region  to  hyperbolic  in  the 
supersonic  region.  In  addition,  special  provision  must  be  made  to 
handle  the  embedded  shock  wave  in  the  flow  field. 

There  are  three  main  categories  of  numerical  methods  for  solving 
steady  inviscid  flow  past  an  airfoil  in  the  transonic  regime.  These 
are  finite  difference  methods,  the  hodograph  method  and  interpolation 
methods . 

Finite  difference  techniques  have  received  the  most  attention  in 
transonic  flow  research  in  recent  years  and  we  now  outline  them 
briefly. 

Magnus  and  Yoshihara  [1970]  first  solved  the  unsteady  Euler 
equations  using  an  explicit  second-order  difference  scheme. 
Unfortunately,  the  method  requires  a  very  large  amount  of  computation 
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time  to  achieve  steady  state  conditions,  and  is  therefore  very  expensive 
for  practical  calculations. 

An  alternative  to  the  time  dependent  approach  is  the  use  of 
relaxation  methods.  Murman  and  Cole  [1971]  successfully  solved  the 
transonic  small  disturbance  equations  by  introducing  a  mixed  finite 
difference  system.  The  direction  of  differencing  is  biased  depending 
upon  whether  the  flow  is  subsonic  or  supersonic.  The  truncation  error 
of  the  difference  scheme  has  the  effect  of  artificial  viscosity,  so 
shock  waves  appear  naturally  during  the  course  of  calculation,  although 
they  are  usually  spread  over  3-4  mesh  points.  The  system  of  difference 
equations  is  solved  by  successive  line  relaxation,  and  the  computed 
results  agree  well  with  experimental  data  for  a  circular  arc  airfoil. 

The  method  was  extended  by  Krupp  and  Murman  [1972]  to  lifting 
airfoils  and  slender  bodies.  Steger  and  Lomax  [1972]  solved  the  full 
potential  equations  for  lifting  airfoils  by  successive  line  over¬ 
relaxation  (SLOR) .  An  interactive  graphic  terminal  is  used  to  change 
the  values  of  circulation  and  relaxation  parameters  as  the  relaxation 
is  proceeding.  To  account  for  flows  not  aligned  to  the  coordinate 
system,  Jameson  [1974]  introduced  a  rotated  differencing  scheme  in 
which  the  direction  of  upwind  differencing  is  rotated  to  conform  with 
the  local  flow  direction.  The  system  of  difference  equations  is  then 
iterated  by  simulating  an  artificial  time  dependent  equation. 

Ballhaus,  Jameson  and  Albert  [1978]  developed  an  implicit  approximate 
factorization  (AF)  algorithm  for  the  solution  of  steady  state  transonic 
small  disturbance  equation,  which  has  a  much  better  rate  of  con¬ 
vergence  than  the  SLOR  algorithm.  Following  the  idea  of  Jameson's 
rotating  upwind  difference  scheme,  Holst  and  Ballhaus  [1979]  solved 
the  full  potential  equation  in  conservation  form  to  ensure  conservative 
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shock  capturing.  Artificial  viscosity  is  added  implicitly  by  retarding 
the  density  according  to  local  Mach  number.  Holst  [1979]  later 
applied  the  method  using  an  arbitrary  mesh,  and  obtained  good  agree¬ 
ment  with  independently  computed  results. 

The  hodograph  method  has  a  long  history  in  transonic  flow  calcula¬ 
tions,  taking  advantage  of  the  property  that  the  governing  equations 
of  plane  motion  become  linear  when  coordinates  in  the  physical  plane 
are  replaced  by  the  velocity  components  as  independent  variables. 

Using  the  hodograph  transformation,  Nieuwland  [1967]  developed  a 
technique  for  computing  shock-free,  symmetrical,  supercritical  flows 
about  quasi-elliptic  airfoil  sections.  The  method  was  later  extended 
by  Boerstoel  [1967]  to  present  a  catalog  of  solutions  for  certain  body 
shapes.  Bauer,  Garabedian  and  Korn  [1972]  also  used  the  hodograph 
method  to  generate  a  shock- free  flow  with  the  corresponding  boundary 
shapes. 

The  final  category  of  the  methods  makes  use  of  quasi-analytic 
techniques.  These  are  the  Method  of  Integral  Relations  (MIR), 

Telenin's  method  and  the  Method  of  Lines  (see  Holt  [1977]).  Each 
approach  uses  smooth  interpolating  functions  to  represent  the  unknown 
variables  in  a  selected  coordinate  direction.  The  partial  differential 
equations  are  thereby  reduced  to  a  set  of  ordinary  differential  equa¬ 
tions  along  a  set  of  rays  in  the  flow  field.  The  resulting  equations 
are  then  solved  as  an  initial  value  problem. 

MIR  was  first  applied  by  Chushkin  [1957]  for  subsonic  critical 
flow  past  an  ellipse  or  ellipsoid.  Later,  Holt  and  Masson  [1971] 
computed  supercritical  flow  about  a  cylinder  with  the  full  potential 
equations.  Tai  [1974]  also  used  MIR  to  solve  the  steady  Euler  equa¬ 
tions  for  a  lifting  airfoil.  Both  of  the  above  methods  located  a 
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shock  point  on  the  body,  but  no  details  about  the  shape  of  the  shock 
in  the  interior  of  the  flow  field  are  obtained. 

Chattot  [1978]  applied  Telenin's  Method  in  the  hodograph  plane 
fbr  flow  past  a  double  wedge.  A  shock  is  fitted  in  the  flow  field  to 
eliminate  the  limit  lines.  The  complete  shock  location  is  obtained, 
but  the  method  is  restricted  to  a  double  wedge,  where  the  boundaries 
in  the  hodograph  plane  are  known  in  advance.  Telenin's  Method  was 
also  used  by  Gross  and  Holt  [1975]  to  calculate  critical  and  super¬ 
critical  shock-free  flow  past  ellipses. 

In  the  present  work  supercritical  flow  past  an  ellipse  at  zero 
angle  of  attack  is  calculated.  The  steady  two-dimensional  full 
potential  equations  are  solved  by  Telenin's  Method  and  the  Method 
of  Lines.  The  jump  conditions  of  the  equations  are  used  to  fit  a 
shock  in  the  flow  field  to  terminate  the  supersonic  region.  The 
formulation  of  the  equations  of  motion  and  the  details  of  the  tran¬ 
sonic  flow  field  are  discussed  in  Section  2.  Applications  of 
Telenin's  Method  and  the  Method  of  Lines  to  the  supercritical  flow 
problem  are  described  in  Section  3.  Section  4  contains  discussions 
of  the  supercritical  calculations.  The  conclusions  are  presented  in 
Section  5. 
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2.0  Formulation  of  the  Problem 

We  consider  the  two-dimensional  flow  of  a  uniform  stream  past  an 
ellipse.  The  free  stream  Mach  number  lies  in  the  high  subsonic  range 
so  that,  while  the  flow  in  the  region  far  from  the  ellipse  is  wholly 
subsonic,  the  flow  field  in  the  neighborhood  of  the  ellipse  is  of 
mixed  type  with  subsonic  regions  near  the  forward  and  rear  stagnation 
points  and  a  local  supersonic  region  near  the  maximum  thickness  section 
of  the  ellipse.  The  local  supersonic  region  is  usually  bounded  by  a 
shock  wave  at  its  downstream  end.  A  typical  flow  pattern  is  shown  in 
Fig.  2(b). 

Viscosity  effects  are  confined  to  the  boundary  layer  near  the 
surface  of  the  ellipse.  The  boundary  layer  calculations  will  be 
carried  out  subsequently  since  these  require  knowledge  of  the  inviscid 
flow  field  as  a  starting  point.  The  shock  terminating  the  local 
supersonic  region  causes  the  boundary  layer  to  separate  so  that  the 
inviscid  and  viscous  flows  interact  significantly.  However,  the  main 
influence  of  boundary  layer  separation  is  to  introduce  an  effective 
thickening  of  the  ellipse  downstream  of  the  shock.  Interaction 
effects  are  therefore  determined  by  integrating  the  inviscid  and 
boundary  layer  equations  separately  and  matching  the  calculations 
along  the  effective  viscous- inviscid  boundary. 

The  shock  wave  introduces  entropy  changes  on  its  downstream 
side.  Provided  that  the  minor-major  axis  ratio  (maximum  thickness 
ratio)  of  the  ellipse  is  sufficiently  small  and,  provided  that  the 
free  stream  Mach  number  is  subsonic,  the  local  Mach  number  ahead  of 
the  shock  wave  will  not  exceed  the  value  1.3.  The  shock  wave  strength 
is  then  sufficiently  small  to  ensure  that  entropy  changes  can  be 
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neglected. 

To  verify  this  statement  we  use  the  following  perfect  gas 
relationship  to  determine  entropy  changes  across  a  shock  wave: 


Y+l 


AS  _  J_  f2yM  -Y+l,  Y  B„  r  (r+l)M 
R  '  Y-l  1 


-]  -  -  tn[- 

Y-l 


-] 


(2.1) 


(Y-l)M  +2 

Here,  M  is  the  incident  Mach  number  (just  upstream  of  the  shock), 

Y  the  ratio  of  specific  heats,  R  the  gas  constant,  and  AS  the  change 
in  entropy.  If  M  =  1 . 3  and  y  =  1.4,  then 


~  =  0.0208 


(2.2) 


Equation  (2.1)  can  be  rewritten  as 

p2 =  Pi(^  )Ytexp(x)]tY’1)  (2-3) 

where  p  is  the  pressure  and  p  the  density.  Subscripts  1  and  2  denote 
conditions  ahead  and  behind  the  shock,  respectively. 

From  Eq.  (2.2) 


[exp(— )](Y'13  =  1.0084  (2.4) 

It  follows  that  for  shocks  with  incidenct  Mach  numbers  M<  1.3  the 
error  introduced  by  the  isentropic  assumption  is  about  0.8?o.  We  can 
therefore  assume,  in  the  present  analysis,  that  the  flow  field  is 
isentropic  and  irrotational . 


2.1  Equations  of  Motion 

We  consider  steady,  two-dimensional  flow  of  a  uniform  air  stream 
past  an  ellipse  at  high  subsonic  free  stream  Mach  numbers.  The  flow 
is  assumed  to  be  inviscid  and  irrotational.  The  governing  equations 


of  motion  are  then 
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Continuity 

Irrotationality 

Writing  these  equations 


The  elliptic  coordinates 
[1972]): 


div(pq)  =  0 

-V 

curl(q)  =  0 

in  elliptic  coordinates: 

3 (puh)  +  9(pvh)  _  Q 
3£  3n 

3(vh)  3 (uh)  =  - 

3£;  *  3n 

C  and  n  are  defined  by 


(2.5) 

(2.6) 

(2.7) 

(2.8) 

ilne-Thomson 


x  =  c  cosh  £,  cos  n  (2.9) 

y  =  c  sinh  £  sin  n  (2.10) 

where  x  and  y  are  Cartesian  coordinates,  c  is  a  constant,  u  and  v 
are  the  velocity  components  in  the  C  and  n  directions,  respectively, 
and  h  is  the  metric  coefficient  given  by: 

h  =  c (sinh^C + sin^n)  2  (2.11) 

It  can  be  shown  that  curves  of  constant  £  represent  confocal 
ellipses  and  curves  of  constant  n  confocal  hyperbolae.  The  foci  of 
the  ellipses  or  hyperbolae  are  located  at  (±c,o).  Some  curves  of 

the  elliptic  coordinate  system  with  c  =  l  are  plotted  in  Fig.  1. 

The  final  equation  to  complete  the  set  is  Bernoulli's  equation: 

H*T02-"o-5\ax2  ‘;i:i 


The  quantity  II  is  the  enthalpy  (per  unit  mass),  q  is  the  flow  speed, 
qmax  bhe  maximum  steady  expansion  speed,  and  the  subscript  o  denotes 
stagnation  conditions.  Equation  (2.12)  may  be  written 


iL«  i-C-iL-)2 

H  lq  > 

o  ^max 


For  a  perfect  gas  with  constant  specific  heats 


(2.13) 


H  = 


-ILL 


(y-i)p 

By  further  assuming  isentropic  flow 


(2.14) 


P  “  P 


(2.15) 


it  follows  that 


JL  =  JL  .  -  f_2_qY-l 

H  p  P  ’ 

o  ‘  o 

Substituting  Eq.  (2.16)  into  (2.13)  and  solving  for  p/pq  we  obtain 


(2.16) 


r-  ■  (i  -  )2i,/v'' 


(2.17) 


•max 


The  three  equations  (2.7),  (2.8)  and  (2.17)  are  the  three  rela¬ 
tions  required  to  determine  the  three  basic  unknowns  u,v  and  p. 

We  now  express  all  variables  in  dimensionless  form  by  dividing 
distances  by  c,  velocities  by  qmax>  and  the  density  by  the  stagnation 
density  p  .  Retaining  the  same  symbols  for  the  non-dimensional 
variables,  the  equations  of  motion  are 


P  = 


3(puh)  3(pvh) 

35  an 

3(vh)  3(uh)  , 

35  “  3n 

(1-q2)1^'1  -  (1-U2-V2)1^-1 


(2.18) 

(2.19) 

(2.20) 
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and 

h  =  (sinh2£  +  sin^n)  1  (2.215 

x  =  cosh  £  cos  n  (2.22) 

y  =  sinh  £  sin  n  (2.23) 

2 . 2  Description  of  the  Flow  Field 

Before  we  proceed  to  solve  the  equations  of  motion,  it  is 
advantageous  to  understand  the  physical  flow  field  and  be  able  to 
choose  an  effective  method  to  solve  the  problem. 

At  very  low  Mach  numbers,  the  compressibility  is  very  small, 
hence  the  flow  can  be  assumed  to  be  incompressible.  The  analytic 
solution  for  incompressible  flow  past  an  elliptic  cylinder  (Milne- 
Thomson  {1972])  is  given  by 

u  =  -jj-  e  cos  n  sinh(£-£Q)  (2.24) 

UO0  S 

v  =  -  -jp  e  0  sin  n  cosh(£-£Q)  (2.25) 

where  is  the  free  stream  velocity  and  £q  the  ellipse  representing 
the  body  of  the  cylinder.  On  the  body,  the  normal  component  u  of 
the  velocity  is  zero,  and  the  tangential  velocity  is  given  by  Eq. 

(2.25).  It  can  be  seen  that  the  flow  accelerates  from  stagnation  at 
the  leading  edge  (n  =  w)  to  a  maximum  speed  at  the  apex  of  the 
cylinder  (n  =  n/2),  and  then  decelerates  back  to  stagnation  at  the 
trailing  edge  (n  =  0).  At  zero  angle  of  attack,  the  flow  is  sym¬ 
metric  about  both  the  y-axis  and  the  x-axis. 

As  the  free  stream  Mach  number  is  increased,  compressibility 
effects  become  more  important,  and  the  compressible  equations  of 


motion  have  to  be  solved.  However,  the  flow  behavior  remains 


qualitatively  the  same  if  the  flow  is  subsonic  throughout  the  field. 
When  the  maximum  local  Mach  number  reaches  the  value  1.0,  the  flow 
is  said  to  be  critical.  The  free  stream  Mach  number  which  produces 
such  a  flow  is  called  the  critical  free  stream  Mach  number. 

For  supercritical  flow,  a  small  supersonic  region  is  embedded 
in  the  subsonic  flow  field.  Although  shock  free  supercritical  flows 
can  be  produced,  they  are  generally  unstable  (see  Busemann  [1949], 
Frankl  [1950],  Guderley  [1953],  Morawetz  [1956,1957])  in  the 
sense  that  a  small  perturbation  of  the  body  contour  in  the  super¬ 
sonic  region  leads  to  a  flow  that  is  discontinuous.  This  will  be 
assumed  the  case  for  our  transonic  flow  calculation.  A  typical 
transonic  flow  field  is  depicted  in  Fig.  2(b). 
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3.0  Numerical  Methods 

Telenin’s  Method  and  the  Method  of  Lines  are  used  as  the 
numerical  schemes  to  solve  the  equations  of  motion.  The  two  methods 
are  very  similar.  In  Telenin's  Method  the  variations  of  the  variables 
in  one  coordinate  direction  are  represented  by  some  smooth  inter¬ 
polating  functions.  In  our  problem,  symmetry  conditions  suggest 
the  use  of  Fourier  series  of  the  form: 


N 

uU ,n)  =  I  a.  (£)cos(i-l)n 
i  =  l  1 


(3.1) 


N-2 

v(£,n)  =  I  b. (C)sin  in 
i  =  l  1 


(3.2) 


where  N  is  the  number  of  rays.  Along  the  jth  ray 

N 

u  =  uU.n.)  =  l  a.  (£)cos(i-l)n. 


J 


i  =  l 


J 


(3.3) 


The  coefficient  a^  can  be  obtained  by  inverting  the  matrix  {cos(i -1 )n^ ) , 

N 

a  =  I  A  u  i  =  1 . N  (3.4) 

j=l  J 

where  {A^ }  =  (cos ( j -1) n^ )  Equation  (3.1)  is  differentiated  to 
obtain  the  n  derivatives,  giving 

8U  au(c.n  )  N 

<!?>t  *  — 3^  ■1£i»i(0(l-i).ln(i-l)n,  (3.5) 

Substituting  Eq.  (3.4)  into  Eq.  (3.5)  yields 


r3u. 


N  N 


%  (1-i)sin(i-1)\ 


(3.6) 


Interchanging  the  order  of  operation, 


,au 


N  N 


N 


*  1  (  E  A.,(l-i)sin(i-l)ne)u.  =  Z^F?.u. 


^  £  j=l  i-1 


(3.7) 


12 


Similarly,  for  the  derivative  of  v 

_  N-2 

•  3v, 


where 


WfiVi  '■'•2 . N 


F.j 

N-2 

G„ .  =  E  B. .  i  cos  in. 

i=l  * 


(3.8) 


(3.9) 


(3.10) 


(B. . }  =  (sin  in .  , ) 

ij  3  +  1 


-1 


(3.11) 


In  the  Method  of  Lines,  the  n  derivatives  are  approximated  by 
finite  differences.  Three-point  or  five-point  difference  schemes  are 
used  depending  on  the  order  of  accuracy  required.  The  derivative 
representation  has  the  same  form  as  given  by  Eq.  (3.7)  with  coef¬ 
ficients  Fjj  derived  from  Taylor  series  expansions.  Hence  in  terms 
of  the  solution  method  and  accuracy,  we  may  consider  the  two  methods 
to  be  equivalent. 

It  is  convenient  to  have  expressions  for  9v/9f;  and  9u/9f;. 

From  Eq.  (2.19) 


Solving  for  9v/9£: 


9h  ,  3v  9h  9u 
v  9?  +  h  H  =  u  9T?  +  h  ^ 


9v  1  ,  9h  .  9u  9h, 
9?  =  h  (U  9^  +  h  9^  '  V 


The  definition  of  h  implies  that 


9h  _  sinh  2E, 
9^  *  2h 


(3.12) 


(3.13) 


and 


(3.14) 
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be  zero,  that  is 


u  =  0  for  £  =  CQ  (3.22) 

where  £q  is  the  elliptic  coordinate  of  the  body.  For  flow  past  an 
elliptic  cylinder  at  zero  angle  of  attack,  the  flow  field  is  sym¬ 
metric  about  the  x-axis,  so 


v  =  0 


and 


^  =  0  for  n  =  0  and  it 
3n 


(3.23) 


Finally,  the  flow  approaches  that  of  a  uniform  free  stream  at  infinity: 


u  -*■  U  cos  n 


and 


v  -+  -  U  sin  n  as  £  -*  °°  (3.24) 

oo 

However,  in  practice,  it  is  more  convenient  to  specify  the  far 
field  boundary  conditions  at  a  finite  distance  from  the  body. 

Following  Murman  and  Cole  [1971],  an  analytical  solution  for  the  far 
field  is  derived  using  transonic  small  disturbance  theory.  The  basic 
transonic  equation  is 


IK*X-  (ylHx/2]x*$_  =  0 

yy 

with  the  variables  and  parameters  defined  by: 


.1/3 

y  *  «  y 


2  ,  2/3 
K  =  (l-M  )/<r* 


(3.25) 


(3.26) 


(3.27) 
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l+62/3q'  + 
U 


(3.28) 


/=  <Sq  7  +  •  • 
U  ^y 

on  * 


(3.29) 


q7  =  <J> 
nx  x 


the  perturbation  velocities 


(3.30) 


Here  6  is  the  thickness  ratio  of  the  airfoil  (or  ellipse),  Mk  the  free 

stream  Mach  number,  q  and  q  are  the  velocity  components  in  the  x  and 

x  y 

y  directions,  respectively. 

We  rewrite  Eq.  (3.25)  in  the  form 


L<t>  =  K*  +*_  =  [  (y  +  1)/2]  (u  )x 

yy 


(3.31) 


Applying  Green's  formula  for  L  in  the  upper  half  plane  and  allowing 
for  a  shock  jump  in  the  flow  field,  we  obtain  the  basic  integral 
equation 


>(x,y)  =  — -  I  - 7 - 2  F(x')dx/ 

it  K  2  J-l  (x-x')"  ♦  Ky 


r+l  1 


~  (x-x7){q'(x',y ')>2 

- ^ dx  7dy 7 


2  2nK*  Ij-ao  (x-x')  +  K(y-y7)2 

The  far  field  is  thus  that  of  the  usual  doublet  for  a  closed  body 


(3.32) 


where 


<f>(x»y) 


2ttK  2  (X  +  Ky  ) 


D  =  doublet  strength 


(5.33) 


-I  i  ( (°°  ? 

=  2  F(x)dx  ♦  {q7(x,y)}Zdxdy  (3.34) 

*  )).<*,  x 

The  doublet  strength  consists  of  the  usual  term  proportional  to 


the  airfoil  volume  and  a  nonlinear  contribution,  unknown  in  advance. 


In  the  numerical  procedure  D  has  to  be  calculated  as  one  of  the 
unknowns  of  the  problem.  Differentiating  Eq.  (3.35), 


/  D  (-x2  ♦  Ky2) 

‘x  r.2  .  ..'2,2 


(3.35) 


2-nK  (x  +  Ky") ' 


/  _  DK2  xy 

V  =  "  u  ,2  ,,-2,2 

7  (x  +  Ky  ) 


(3 . 36) 


The  flow  velocities  expressed  in  the  £  and  n  directions  are  given  by: 


u  =  jj-  (qx  sinh  4  cos  n  +  q  cosh  £  sin  n) 


(3.37) 


v  =  ^-  (-qx  cosh  £  sin  n  +  qy  sinh  £  cos  n) 


(3.38) 


Substitution  of  Eqs.  (3.28)  and  (3.29)  into  (3.37)  and  (3.38)  gives 
the  necessary  boundary  conditions  in  the  far  field. 


3.2  Jump  Conditions 

For  transonic  flow  at  sufficiently  high  subsonic  free  stream 
Mach  number,  the  flow  becomes  supercritical.  A  region  of  local 
supersonic  flow  is  developed  over  the  maximum  thickness  region  of  the 
body  and  this  is  terminated  on  the  downstream  side  by  a  shock  wave. 

In  the  inviscid  flow  approximation,  the  shock  wave  is  modeled  by  a 
jump  discontinuity  in  the  solution.  To  ensure  uniqueness,  we  require 
that  entropy  increases  across  the  shock  wave.  For  the  full  potential 
approximation,  the  entropy  change  is  assumed  to  be  negligible. 
Uniqueness  is  attained  by  allowing  only  the  existence  of  compression 
shocks,  but  not  expansion  shocks.  The  jump  conditions  for  the  full 
potential  equations  are  different  from  the  usual  Rankine-Hugoniot 
relations,  and  can  be  derived  by  writing  the  equations  of  motion  in 
conservation  form.  Applying  the  two-dimensional  form  of  the  divergence 
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theorem  to  Eqs.  (2.18)  and  (2.19),  we  obtain 


and 


<puh>(dn)s  -  <pvh>(d£)s  =  0 


(5.39) 


<vh>(dn)s +  <uh>(d£)s  =  0  (3.40) 

where  <  >  denote  a  jump  in  the  quantity  across  the  shock  and  subscript 
s  denotes  an  element  in  the  shock  surface.  Equations  (3.39)  and  (3.40) 
can  be  rewritten  as 


where 


<puh>n^  -  <pvh>  =  0 

<vh>n  '  +  <uh>  =  0 
s 


ns  =  (^s  *s  t,ie  s*10ck  wave  ar>gle 


(3.41) 

(3.42) 

(3.43) 


In  Eqs.  (3.41)  and  (3.42)  h  is  the  metric  coefficient  which  depends 
only  on  the  geometry  of  the  coordinate  system  and  is  continuous 
throughout  the  field,  hence  can  be  eliminated  from  Eqs.  (3.41)  and 
(3.42).  The  final  form  of  the  jump  conditions  is  then 


<pu>Hg  -  <pv>  =  0  (5.44) 

<v>n^  +  <u>  =  0  (3.45) 

Equations  (3.44)  and  (3.45)  represent,  respectively,  the  conservation 
of  mass  flux  and  continuity  of  tangential  velocity  across  the  shock 
wave.  The  density  p  in  Eq.  (3.44)  is  given  by  Bernoulli's  equation 
(2.20).  Thus  the  jump  conditions  for  the  full  potential  equations  arc 
completely  specified. 
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It  is  interesting  to  compare  the  shock-fitting  and  shock¬ 
capturing  methods.  In  the  finite  difference  treatment  of  super¬ 
critical  flow,  artificial  viscosity  is  added  to  the  differential 
equations  as  a  result  of  the  truncation  errors  generated  by  the 
difference  equations.  No  explicit  jump  conditions  are  needed, 
provided  that  the  equations  are  written  in  divergence  form  to 
conserve  mass  flux.  Shock  waves  evolve  naturally  during  the  course 
of  the  calculation,  although  they  usually  spread  over  several  mesh 
points.  In  principle,  the  shock  wave  can  be  made  arbitrarily 
sharp  by  refining  the  mesh  points  near  to  it;  however,  this  will 
slow  down  the  rate  of  convergence  considerably.  By  employing  a 
shock-fitting  technique,  the  jump  conditions  are  satisfied  exactly. 

The  shock  wave  is  perfectly  sharp,  hence  no  refinement  is  necessary. 
The  drawback  of  this  method  is  that  the  iteration  may  not  converge 
if  the  initial  guess  of  the  shock  location  is  too  inaccurate. 

3.3  Singular  Points 

In  Section  3.0  we  assumed  a  Fourier  series  representation  in  the 
n  direction  for  the  unknown  flow  quantities,  and,  as  a  result, 
derived  a  set  of  first  order  ordinary  differential  equations  (Eqs. 
(3.16)  and  (3.19))  in  the  E;  direction.  However,  it  is  also  pos¬ 
sible  to  assume  an  analytic  representation  in  the  e;  direction  and 
obtain  a  set  of  ordinary  differential  equations  in  the  n  direction. 

The  advantages  and  disadvantages  of  each  formulation  will  become 
apparent  at  a  later  stage. 

For  the  latter  formulation,  we  must  derive  expressions  for 
3u/9q  and  9v/9n.  Rearranging  Eqs.  (2.18)  and  (2.19),  and  solving 
for  3v/9n  and  9u/9n,  we  obtain 


1.) 


3u 

3n 


3v  v  sinh  2K  u  sin  2n 

3C  +  ^2  '  .,2 


2h 


2h 


(3.46) 


3v  _  ^2_ 

3n  '  Q, 


(3.47) 


where 


n  _.  ,3u  dv  .  2  3u 

P2  =  ^uv(— +  — )  +  2u  w 


,  ....  2  2W8U  u  sinh  2£  v  sin  2ns  ,, 

-  (y-1)  (1-u  -v  )  (gjr  +  - 5 - -  +  - - - )  (3.48) 


2h 


2h‘ 


Q,  =  (y-1)(1-U2)  -  (Y+l)v2 


(3.49) 


When  Eqs.  (3.16)  and  (3.19)  are  examined  in  detail,  it  is 
observed  that  they  have  a  saddle  point  singularity  when  the  denominator 
Qj  becomes  zero,  that  is, 

(Y-l) (1-v2)  -  (y  +  1)u2  =  0  (3.50) 


after  rearranging,  we  obtain 


Y  +  l 


2 


v 


1 


(3.51) 


which  represents  an  ellipse  in  the  u,v  plane.  From  Bernoulli's 
equation: 


q2  =  (y-DM2 

2  +  (y-1)M2 

* 

so  the  non-dimensional  critical  velocity  q  is 


(3.52) 


"2  _  rl 

"  Y+1 


(3.53) 


Substituting  Eq.  (3.53)  into  (3.51),  the  ellipse  of  singularities  may 


be  written  as 


20 


2 

u 


+ 


1 


(5.54) 


The  singular  ellipse  and  the  sonic  circle  are  both  plotted  in  Fig.  3a. 

We  can  see  that  all  points  on  the  ellipse  lie  outside  the  sonic  circle, 

* 

except  for  v = 0,  u =  ±q  .  For  critical  flow,  only  one  point  is  on  the 

* 

sonic  circle,  namely  u  =  0,  v  =  ±q  ,  therefore  there  are  no  singularities 
for  a  critical  flow  calculation  when  integrating  away  from  the  body. 

It  is  apparent  that  no  singularities  will  be  encountered  even  for 
supercritical  flow  calculations. 

On  the  other  hand,  using  the  second  formulation  and  integrating 
in  the  n  direction,  th .  denominator  of  Eq .  (3.47)  becomes  zero  when 


(Y-l)(l-u2)  -  (y+1)v2  =  0 


(3.55) 


Hence  the  singular  ellipse  in  the  u,v  plane  is  given  by 


u 


2 


1 


(3.5b) 


which  is  shown  in  Fig.  3b.  But  any  critical  or  supercritical  flow  has 

* 

a  point  on  the  body  with  u  =  0  and  v  =  -q  ,  which  is  a  point  on  the 
ellipse  given  by  (3.56).  Therefore  it  is  obvious  that  integration 
in  the  n  direction  always  leads  to  at  least  one  singularity  at  sonic 
points  on  or  near  the  body. 


3.4  Implementation  of  the  Numerical  Scheme 

As  discussed  in  the  preceding  section,  saddle  point  singularities 
will  arise  if  we  assume  interpolating  functions  in  the  £;  direction. 
Hence  it  will  be  appropriate  to  use  the  first  formulation  given  in 
Section  3.0.  Expressions  (3.7)  and  (3.8)  are  substituted  into  Eqs . 
(3.16)  and  (3.19)  to  form  a  set  of  (2*N  -  2)  simultaneous  ordinary 
differential  equations.  At  £  =  the  flow  tangency  condition  on 


the  body  is  given  by  Eq.  (3.22).  An  initial  estimate  of  values  of 
the  tangential  velocities  vq  on  the  surface  is  made,  and  using  these 
as  initial  data,  the  equations  are  integrated  away  from  the  ellipse 
£  =  CQ.  A  variable  step,  fifth-order  Runge-Kutta  method  is  used  to 
integrate  the  differential  equations.  The  integration  is  terminated 
at  a  distance  sufficiently  far  away  from  the  body,  say  £ =  £  ,  which 
will  be  defined  later.  The  velocities  v  calculated  at  £  are  then 

oo 

compared  with  the  far  field  velocities  given  by  Eq.  (3.38).  If  the 
two  sets  of  values  differ,  tangential  velocities  on  the  surface  are 
then  adjusted  and  the  integration  is  repeated.  The  procedure  is 
repeated  until  the  far  field  solution  converges.  This  can  be  done 
very  efficiently  by  the  use  of  Powell's  method  (Powell,  1964),  which 
minimizes  the  sum  of  squares  of  the  differences  between  the  far  field 
velocities  by  adjusting  the  surface  velocities.  The  iteration  is 
terminated  when  the  sum  satisfies  the  specified  tolerance. 

Gilinskii,  Telenin  and  Tinyakov  [1964]  pointed  out  that  solving 
a  Dirichlet  problem  as  a  Cauchy  problem  is  inherently  unstable  with 
respect  to  the  prescribed  data.  This  phenomenon  is  known  as 
Hadamard  instability.  Jones  and  South  [1972]  also  encountered 
Hadamard  instability  in  applying  the  Method  of  Lines  and  found  growth 
in  error  proportional  to  exp(N£),  where  N  is  the  number  of  rays  and  £ 
the  direction  of  integration.  As  a  consequence,  the  number  of  rays 
and  the  far  field  distance  have  to  be  restricted.  However,  if  we  wish 
to  obtain  a  solution  with  reasonable  accuracy,  we  must  employ  suf¬ 
ficient  number  of  rays  to  represent  the  variables.  The  occurrence 
of  a  shock  wave  on  the  body  makes  it  even  more  desirable  to  have  a  de¬ 
tailed  representation  near  the  body. 

To  overcome  the  above  difficulties,  we  propose  to  solve  the 
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problem  in  two  stages.  In  the  first  stage,  a  very  coarse  representation  j 

i 

of  the  variables  is  used,  which  enables  us  to  integrate  the  equations  I  ' 

away  from  the  body  to  the  far  field  without  instability  problems.  A 
supercritical  shock  free  flow  is  obtained  from  this  calculation. 

However,  as  discussed  in  Section  2.2,  shockless  flows  are  known  to  be 

unstable  and  not  likely  to  occur  in  practical  situations.  Hence,  in  j 

our  supercritical  flow  calculation,  we  always  assume  that  the  flow  is  j- 

discontinuous.  In  order  to  model  the  shock  wave,  we  have  to  treat  the  j- 

region  near  the  body  in  a  different  manner.  The  procedure  for  the 
different  stages  is  described  in  the  following  sections. 

3.4.1  Coarse  Solution  • 

In  this  stage,  we  assume  Fourier  series  representation  of  the  [ 

form  (3.1)  and  (3.2)  for  u  and  v.  For  flow  over  an  ellipse,  we  further 

notice  that  the  flow  is  symmetric  about  the  y-axis  when  the  solution  1 

is  smooth.  Thus  we  can  economize  on  the  number  of  rays  by  assuming 
the  series  of  the  form: 

N-l  I 

u(f.,n)  =  l  a.  (f.)cos(2i-l)n  (3.57) 

i  =  l  1  i 

N-l  j 

v(£,n)  =  Z  b. (C)sin(2i-l)n  (3.58) 

i  =  l  1  j 

So  for  the  coarse  solution,  we  only  need  to  compute  the  flow  in  the 

i 

second  quadrant.  Equations  (3.16)  and  (3.19)  are  integrated  [ 

simultaneously  from  CQ  to  C^.  Boundary  conditions  are  satisifed  by 
using  the  procedure  described  in  Section  3.4.  At  supercritical  Mach 

number,  a  continuous  flow  with  a  small  embedded  supersonic  zone  is  j 

obtained,  and  is  depicted  in  Fig.  2a.  j 
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3.4.2  Refined  Solution  Near  the  Body 


Although  the  coarse  solution  does  not  have  enough  accuracy  to 
resolve  the  shock  wave  which  occurs  near  the  body,  it  provides  a 
fairly  good  representation  of  the  flow  field  away  from  the  body 
where  the  flow  is  smooth.  The  strategy  here  is  to  use  a  larger 
number  of  rays  to  represent  the  flow  field  close  to  the  body,  the 
coarse  solution  at  an  intermediate  value  of  £,  say  is  used  as 
the  outer  boundary  condition  for  the  refined  solution  near  the  body. 

In  this  way,  the  distance  in  the  £  direction  is  kept  small  and  a 

larger  number  of  rays  can  be  used  without  causing  instabilities. 

It  is  advantageous  to  integrate  the  equations  in  the  n  direction 
if  we  wish  to  fit  a  shock  in  the  flow  field.  Following  the  idea  of 
Fletcher  [1975],  we  divide  the  region  near  the  body  into  two  parts. 
The  forward  part  is  enclosed  by  C  =  S0»  £  =  ^  ,  n  =  u,  and  n  =  n^,  and 

the  rear  portion  by  £  =  4q,  £  =  £^ ,  n  =  and  n  =  0,  which  we  shall  call 

region  1  and  region  2,  respectively.  The  configurations  are  shown 
schematically  in  Figs.  4  and  5.  f  is  chosen  such  that  the  point 
(£^,n^)  is  at  the  top  of  the  sonic  line  (see  Fig.  2a).  In  our  case, 
ni  =  n/2. 

In  region  1,  finite  difference  formulae  of  the  form  (3.7)  and 

(3.8)  are  used  to  calculate  the  derivatives  in  the  n  direction.  The 

ordinary  differential  equations  (3.16)  and  (3.19)  are  integrated  as 

in  the  coarse  calculation.  No  boundary  condition  is  required  on 

n  =  tt/2,  since  there  is  no  influence  from  downstream  in  a  supersonic 

region,  smooth  transition  through  the  sonic  line  will  be  sufficient 

for  uniqueness.  Tangential  velocities  on  the  body  are  adjusted  so 

that  velocities  v  at  £.  match  those  of  the  coarse  calculation.  Five 

i 

rays  have  been  used  in  this  region  without  encountering  stability 
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problems . 

In  region  2,  derivatives  in  the  £  direction  are  calculated  by 
finite  difference  formulae  of  the  following  form: 


.  M 

=  1  E»-u- 
'■3£ J  l  .  ,  i.1  l 
1=1 


„  M 

,3  V.  _  r  u 

3C  £  "i=1H£iVi 


(3.59) 


(3.60) 


where  M  is  the  number  of  rays,  E^  and  are  matrices  obtained  by 
Taylor  series  expansion.  Expressions  (3.59)  and  (3.60)  are  substituted 
into  Eqs.  (3.46)  and  (3.47)  to  form  a  set  of  ordinary  differential 
equations.  At  r)  =  v/2,  the  converged  solution  from  region  1  is  used  as 
the  initial  condition  to  integrate  the  equations  from  n = u/2  to  n = 0. 

At  n = 0,  the  symmetry  condition  v = 0  is  imposed.  However,  in  the 
supersonic  region  the  flow  has  no  forewarning  of  the  downstream  con¬ 
ditions  and  the  flow  will  not  be  able  to  adjust  to  satisfy  symmetry 
conditions  at  n=0.  Physically,  the  supersonic  region  is  terminated  by 
a  shock  wave,  the  subsonic  region  behind  the  shock  wave  is  subsequently 
compressed  to  satisfy  the  boundary  condition  downstream.  To  account 
for  the  embedded  shock  wave,  the  equations  are  integrated  to  an 
intermediate  value  of  n,  say  ns(C)»  across  which  jump  conditions 
(3.44)  and  (3.45)  are  applied.  The  integration  is  then  resumed  and 
carried  out  until  n = 0  is  reached.  Powell's  method  is  used  to  adjust 
the  shock  location  ng  until  v  becomes  zero  on  n  =  0.  The  location  of 
the  shock  is  specified  by  its  location  on  the  surface  and  the  shock 
slope  on  other  rays.  It  is  known,  a  priori,  that  the  local  shock 
shape  must  be  normal  to  the  surface  in  order  to  preserve  the  boundary 
condition  of  zero  normal  flow  at  the  surface. 

By  splitting  the  solution  domain  near  the  body  into  two  parts,  we 
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have  been  able  to  integrate  the  equations  of  motion  in  different 
directions.  In  the  rear  part  we  chose  to  integrate  the  equations  in 
the  n  direction  so  that  a  shock  can  be  fitted  in  the  flow  field. 

Since  the  flow  is  supersonic  ahead  of  the  shock  and  subsonic  behind 
(at  least  when  boundary  layer  interaction  effects  are  not  considered) , 
no  saddle  point  singularity  is  encountered,  thus  the  integration  can 
be  carried  out  without  difficulty.  As  a  word  of  caution,  we  note 
that  the  shock  wave  does  not  extend  all  the  way  from  the  surface  to 
£,  =  (see  Fig.  2b).  If  we  use  a  large  number  of  rays,  the  rays  of 
constant  £;  with  values  close  to  ^  will  pass  through  the  sonic  line 
and  will  cause  difficulty  if  we  try  to  integrate  through  this  line. 
Therefore  care  must  be  taken  to  ensure  that  no  rays  pass  through  the 
sonic  line. 

After  a  converged  solution  for  region  2  is  obtained,  the  doublet 
strength  D  will  be  re-calculated,  and  the  whole  procedure  repeated. 

The  solution  is  considered  to  have  converged  globally  when  the  values 
of  D  at  successive  iterations  agree  to  within  the  prescribed  tolerance. 
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3.4.3  Powell's  Method  ! 

I 

To  complete  the  description  of  our  numerical  scheme,  we  shall  | 

describe  Powell's  method  briefly;  a  more  detailed  analysis  can  be  found 
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M  N  3e  3e.  N  3e. 

,l  {  1  3F_3F7}6Fj  =  ~  .Z  Ck  ‘ST 
j=l  k=l  1  j  J  k=l  i 


(3.62) 


New  values  of  F  are  given  by 


F  =  F  ,  ,  +  A  6F 
-  -old  — 


(3.63) 


in  which  A  is  chosen  (by  search)  such  that  is  minimized  along  the 

direction  6F.  During  the  search  along  ^F  to  locate  the  minimum, 

functions  e.  have  to  be  evaluated  at  different  values  of  A:  thus  one 
i 

can  calculate  the  rate  of  change  of  along  the  direction  at  the 
new  minimum  point  by  finite  differences.  Powell  shows  how  these 

partial  derivatives  can  be  used  in  conjunction  with  previous  values 

3ck 

of  jpr-  to  determine  values  for  the  next  step  given  by  (3.62). 
i 

In  principle  the  method  guarantees  convergence  since  a  step  is 
2 

taken  only  when  Ze^  decreases.  It  also  has  quadratic  convergence 
provided  one  is  sufficiently  near  the  solution  and  =  0  at  the 


minimum. 


4.0  Results  and  Discussion 


The  algorithm  introduced  in  the  previous  section  is  evaluated 
in  this  section  by  presenting  a  range  of  numerically  computed  solutions. 
Ellipses  with  thickness  ratio  6  =  0.4  and  6  =  0.2  are  chosen  for  the 
test  cases.  Free-stream  Mach  numbers  are  assumed  to  be  high  enough 
so  that  a  shock  wave  will  always  occur.  Gross  and  Holt  [1975] 
reported  critical  flow  for  6  =  0.4  at  M^  =  0.587.  Symmetric,  super¬ 
critical,  shock-free  flows  were  obtained  up  to  M  = 0.644.  A  range  of 
free-stream  conditions  has  been  chosen  for  our  computations.  For 
thickness  ratio  6  =  0.4,  free-stream  Mach  numbers  were  chosen  to  be 
0.65,  0.66,  0.67  and  0.68.  For  6  =  0.2,  M  =0.77,  £  is  assumed  to 

OO  CO 

be  2.5,  which  was  found  to  be  sufficiently  large  by  Gross  and  Holt 

[1975].  Three  rays  are  used  for  the  coarse  calculation.  When  four 

rays  are  used,  the  solution  tends  to  oscillate  in  the  £  direction  at 

large  £ ,  which  is  due  to  the  instability  discussed  in  Section  3.4. 

For  region  1  of  the  refined  calculation,  five  rays  are  used  without 

encountering  instability  problems.  However,  for  region  2  of  the 

refined  calculation,  only  three  rays  can  be  used.  Besides  the 

instability  problem,  there  are  possible  singular  points  in  the 

ordinary  differential  equations  depending  on  whether  or  not  the  rays 

pass  through  the  sonic  points.  This  is  due  to  the  fact  that  the  shock 

does  not  extend  all  the  way  out  to  £  =  f^,  typically,  the  top  of  the 

shock  is  located  at  approximately  two-thirds  of  the  distance  between 

£  and  £..  It  follows  that  the  use  of  more  than  three  rays  will  cause 
o  1 

at  least  one  ray  to  pass  through  the  sonic  point,  which  is  hazardous 
when  integrating  in  the  o  direction. 

The  present  results  are  compared  with  calculations  using  the 


28 


shock-capturing  method  of  Holst  [1979]  and  are  shown  in  Figs.  6-10. 

For  6  =0.4,  at  M  =0.65,  0.66  and  0.67,  the  two  methods  agree  very- 
well  with  the  shock  locations  on  the  body  almost  identical.  The 
surface  velocity  profiles  obtained  by  both  methods  show  similar 
characteristics.  The  flow  undergoes  a  small  compression  before  the 
shock  wave  is  encountered.  Behind  the  shock,  a  small  post  shock 
expansion  wave  is  observed,  after  which  the  flow  is  recompressed  back 
to  stagnation  condition  at  the  trailing  edge.  The  two  methods  show 
the  largest  discrepancies  near  the  shock  wave;  in  all  the  three  cases 
tested,  Holst's  method  consistently  obtains  a  higher  maximum  velocity- 
on  the  body  and  shows  a  steeper  pre- shock  compression.  At  M  =0.64, 
a  solution  with  an  embedded  shock  wave  could  not  be  obtained  by  the 
present  method.  As  can  be  seen  in  the  solution  at  Ma  =  0.65,  the  shock 
jump  is  very  weak,  and  it  is  quite  probable  that  symmetric  shock-free 
flows  exist  for  free-stream  Mach  numbers  lower  than  0.65.  However, 
using  Holst's  method,  a  solution  with  shock  jump  is  obtained  for 
M  =0.64.  At  M  =0.68,  the  solution  shows  a  local  Mach  number  of 

OO  00 

1.61  ahead  of  the  shock,  the  isentropic  assumption  at  this  Mach  number 
will  introduce  an  error  of  about  4.7%  according  to  Eqs.  (2.3)  and  (2.5). 
Hence  any  solutions  obtained  at  or  above  this  free-stream  Mach  number 
will  be  erroneous. 

Due  to  the  unstable  nature  of  the  present  method,  the  round-off 
error  grows  as  exp(N*£),  where  N  is  the  number  of  rays  in  the  £,  direc¬ 
tion.  It  is  difficult  to  assess  the  accuracy  of  the  method.  Neverthe¬ 
less,  good  agreement  is  obtained  between  the  present  method  and  the 
shock-capturing  method.  The  present  method  takes  about  3  seconds  to 
execute  on  a  CDC  7600,  whereas  Holst's  method  takes  about  6  seconds  for 
a  mesh  size  of  90*40. 


Unfortunately,  the  present  method  does  not  guarantee  convergence 
unless  the  initial  guess  is  reasonably  close  to  the  converged  solution 
One  remedy  is  to  increase  the  free- stream  Mach  number  by  a  small  frac¬ 
tion  at  a  time,  say  by  0.005,  and  then  use  the  solution  obtained  for 
a  previous  Mach  number  as  an  initial  guess.  Otherwise,  good  judgment 
and  trial  and  error  are  needed  in  providing  a  good  initial  guess. 
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5.0  Conclusions 


A  composite  numerical  scheme  has  been  developed  which  is  based 
in  part  on  Telenin's  Method  and  in  part  on  the  Method  of  Lines.  The 
numerical  method  has  been  designed  to  solve  for  supercritical  flow 
over  an  ellipse,  when  the  free-stream  Mach  number  is  high  enough  to 
generate  an  embedded  shock  wave  in  the  flow  field.  A  fitting 
technique  is  used  to  determine  this  shock  so  that  the  Rankine- 
Hugoniot  jump  conditions  are  satisfied  exactly  across  the  shock  wave. 

Good  agreement  is  obtained  between  the  present  method  and  the 
shock-capturing  technique.  Further  improvement  in  the  present  method 
can  be  achieved  by  the  introduction  of  non-svmmetrical  flow  effects 
into  the  solution  over  the  forward  part  of  the  ellipse.  To  this  end, 
it  is  desirable  to  represent  the  far  flow  field  solution  in  terms  of 
distributed  singularities  along  the  major  axis  of  the  ellipse  rather 
than  in  terms  of  singularities  all  located  at  the  ellipse  center. 
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FIG.  2(a)  UNSTABLE  SUPERCRITICAL  SHOCK  FREE  FLOW 
FIELD.  UPPER  HALF  PLANE. 


FIG.  2(b)  TYPICAL  SUPERCRITICAL  FLOW  FIELD.  UPPER 
HALF  PLANE. 
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f>.  3(a)  SINGULAR  ELLIPSE  FOR  FIRST  FORMULATION. 
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FIG.  3(b)  SINGULAR  ELLIPSE  FOR  SECOND  FORMULATION, 
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BOUNDARY  OF  LOCAL  SUPERSONIC  REGION.  6=0.4,  M  =0.66 


